CompGen2021: The Capstone Project

Quick Start

All data can be accessed at here.

• You can download the data by using Linux/Ubuntu command window as well. To download the data, copy the link of the file from the HTML file provided to you and just type wget link.

• For example: wget https://bimsbstatic.mdc-berlin.de/akalin/AAkalin_DrugResponse/prepared/AAkalin_DrugResponse-prepared.tar.gz. Typing this command will download the entire data folder.If you’re having any issues on wget command, you can get more information from the link here.

• Once you download the data, you need to unzip it.You can use gunzip() function to unzip the file with .gz extension. Additionally, you may also need to use R.utils::untar() function to untar the .tar files. All of the data is in .tsv format. There is no significantly difference between the .csv and .tsv format, so you can simply read .tsv by specifying sep = '\t', or by using read.delim() function.

• For example: read.delim(R.utils::gunzip("data/CCLE/drug_response/TRAMETINIB.tsv.gz")) and you can assign the data you read to a variable.

From the code below, you can see how the PCA graph of gene expression is plotted by coloring it based on the drug.

#read the gene exp data
CCLE_gex <- read.delim("data/CCLE/all_intersect_genes/gex.tsv") %>% 
  rename(cell_lines = X) #Change the name of cell lines column to cell_lines
# with the dimension of 1305 x 18356

#read the drug response data
drug_resp_CCLE <- read.delim("data/CCLE/drug_response/ALPELISIB.tsv") %>%
    dplyr::bind_rows(read.delim("data/CCLE/drug_response/BUPARLISIB.tsv")) %>%
    dplyr::bind_rows(read.delim("data/CCLE/drug_response/CETUXIMAB.tsv")) %>%
    dplyr::bind_rows(read.delim("data/CCLE/drug_response/ERLOTINIB.tsv")) %>%
    dplyr::bind_rows(read.delim("data/CCLE/drug_response/RIBOCICLIB.tsv")) %>%
    dplyr::bind_rows(read.delim("data/CCLE/drug_response/TRAMETINIB.tsv")) %>%
    rename(drug = column_name, cell_lines = sample_id)
#with the dimension of 4190 x 5

#Remove near zero variation 
nzv <- caret::preProcess(CCLE_gex,method="nzv",uniqueCut = 15)
CCLE_gex_pp <- predict(nzv,CCLE_gex)
#with the dimension of 3045 x 16023

#center scale and pca
preprocessParams  <- caret::preProcess(CCLE_gex_pp,method=c("pca"))
CCLE_gex_pp <- predict(preprocessParams,CCLE_gex_pp)

#combine the drug response values
ccle_data <- dplyr::inner_join(drug_resp_CCLE, CCLE_gex_pp, by = "cell_lines") %>%
  select(!data_type  & !variable_type)
#with the dimension of 3045 x 18358

#Save the plots
ALPpcaplot <- ccle_data %>% filter(drug == "ALPELISIB") %>%
  ggplot(aes(x=PC1,y=PC2))+
  geom_point(aes(color =value))+
  ggtitle("ALPELISIB")+
  theme(legend.position = "right")+
  scale_color_gradient(low = "orange", high = "Blue")

BUPpcaplot <- ccle_data %>% filter(drug == "BUPARLISIB") %>%
  ggplot(aes(x=PC1,y=PC2))+
  geom_point(aes(color =value))+
  ggtitle("BUPARLISIB")+
  theme(legend.position = "right")+
  scale_color_gradient(low = "orange", high = "Blue")
    
CETpcaplot <- ccle_data %>% filter(drug == "CETUXIMAB") %>%
  ggplot(aes(x=PC1,y=PC2))+
  geom_point(aes(color =value))+
  ggtitle("CETUXIMAB")+
  theme(legend.position = "right")+
  scale_color_gradient(low = "orange", high = "Blue")

ERLpcaplot <- ccle_data %>% filter(drug == "ERLOTINIB") %>%
  ggplot(aes(x=PC1,y=PC2))+
  geom_point(aes(color =value))+
  ggtitle("ERLOTINIB")+
  theme(legend.position = "right")+
  scale_color_gradient(low = "orange", high = "Blue")

RIBpcaplot <- ccle_data %>% filter(drug == "RIBOCICLIB") %>%
  ggplot(aes(x=PC1,y=PC2))+
  geom_point(aes(color =value))+
  ggtitle("RIBOCICLIB")+
  theme(legend.position = "right")+
  scale_color_gradient(low = "orange", high = "Blue")

TRApcaplot <- ccle_data %>% filter(drug == "TRAMETINIB") %>%
  ggplot(aes(x=PC1,y=PC2))+
  geom_point(aes(color =value))+
  ggtitle("TRAMETINIB")+
  theme(legend.position = "right")+
  scale_color_gradient(low = "orange", high = "Blue")

cowplot::plot_grid(ALPpcaplot,BUPpcaplot,CETpcaplot,ERLpcaplot,RIBpcaplot,TRApcaplot)

• NOTE : To use mutation data, you need to make it a binary matrix. However, you may have limited memory problems with using this function. We recommend you run the code and save the data as .rds, then read the .rds file in a new session. You can do it by following the code snippet below.

#necessary libraries
library(reshape2)
library(data.table)

#read the data
CCLE_mut_pan <- read.delim("data/CCLE/oncoKB_intersect_genes/mut.tsv") %>%
  select(sample_id,gene_id,protein_change) %>% #Select the necessary columns
  as.data.table() %>% melt(id.vars = c("sample_id","gene_id"))

CCLE_mut_pan <- dcast(CCLE_mut_pan, sample_id ~ gene_id + value, length)
saveRDS(CCLE_mut_pan,"CCLE_mut_pan.rds") #save the data as a RDS file

#Now open a new session then you can read and assign the data as a data frame
ccle_mut_pan <- as.data.frame(readRDS("ccle_mut_pan.rds"))

Get to Know The Data

knitr::include_graphics(c("PDX.png","cellculture.jpg"))

TASK 1: Use Multi-Omics Data to Predict Drug Response

This task is about using multi-omics data to predict drug response (gene expression, CNVs and mutations). The datasets to be used are in all_intersect_genes and drug_response folders in GDSC_Broad, GDSC_Sanger, CCLE and PDX_all folders. PDX_all should not be used for training, only for testing. The model parameters during training should not be optimized in anyway towards PDX prediction performance. This is extremely important. You should use all omics dataset in training not just expression.

The sub-tasks are as follows:

0. Check if you need some sort of normalization, filtering and/or transformation. These methods have to be applied in all training and test datasets. More complicated methods for batch effect correction is on Task 3 and don’t belong in this task. Just check if normalization, transformation and filtering is needed. And apply whatever you think is necessary. You should first look at normalization and filtering per omic data type. You then can choose to concatanate filtered and normalized single omic datasets into a multi-omics matrix prior to training. You can also experiment to see how different data processing techniques affects accuracy. But those experiments do not belong in the final report, just include the final processing method and state if you did some experimentation to arrive at that method.

1. Train machine learning models to predict drug response on each cell line dataset (GDSC_Sanger, GDSC_Broad, CCLE) using 10-fold cross-validation. You can use any machine learning model available to you in R and tune its performance. You would train a model for each drug and each cell line dataset. For those best, performaing models, make a table or heatmap (with values printed on the cells) for the CV R-squared values. columns: CCLE, GDSC_Broad, GDSC_Sanger and rows: ALPELISIB, BUPARLISIB, CETUXIMAB, ERLOTINIB, RIBOCICLIB, TRAMETINIB. Do not change this row/column order.

2. By now you have a tuned model for each drug and cell line dataset combination. What happens when you use a model trained on one cell line dataset on on the other two? For example, for a given drug say Cetuximab, apply the CCLE model on GDSC_Sanger and GDSC_Broad datasets and record R-squared for those predictions. Then, you can take the GDSC_Broad model and apply it on other two and record R-squared values. If you keep doing this, for each drug you will have six values for the R-squared measure. Output median and mean values as a table, and make boxplots for R-squared values of each drug as obtained above. Is there a dataset that results in poor models for many drugs?

3. Take the best model from sub-task 1 and also the combine the models models in an ensemble learner ( this could be simply averaging their output) and predict drug response on the PDX data set. Plot the metrics (R-squared, Pearson correlation and Spearman’s Correlation) for observed vs predicted efficacy of each drug. Did the simple model perform better than the ensemble model?

• HINT : How to decide the best model from step-1 for a given drug, this could be simply looking at CV R-squared values, or pick the model which has similar CV R-squared values to the out-of-sample R-squared values on the two other datasets. Anything reasonable is OK.

• HINT : For the ensemble model, you can simply average output or re-train an ensembl model with cross-validation. You don’t have to use all the models in the ensemble, you can choose to remove a model if you believe it performs poorly. Anything is OK but not every approach will be performant. Use your imagination.

• IMPORTANT NOTE : the PDX data set is for testing ONLY and it shouldn’t be used for training the model or optimizing any training parameter including ensemble models.

TASK 2: Predict Drug Response Using Panel Sequencing Data.

This task is about using panel sequencing data to predict drug response (CNVs and mutations for few hundred genes). This kind of data is usually obtained as part of cancer diagnostics therefore more easily accesible for clinical samples. However, since the panel sequencing data measures limted set of features from a limited set of genes, it may not be the best choice for clinical variable modeling. The datasets should be used for this task are in oncoKB_intersect_genes and drug_response folders in GDSC_Broad, GDSC_Sanger, CCLE and PDX_all folders. PDX_all should not be used for training, only for test. The model parameters during training should not be optimized in anyway towards PDX prediction performance. This is extremely important.

The sub-tasks are as follows:

0. Check if you need some sort of normalization, filtering and/or transformation. These methods have to be applied in all training and test datasets. More complicated methods for batch effect correction is on Task 3 and don’t belong in this task. Just check if normalization, transformation and filtering is needed. And apply whatever you think is necessary. You can also experiment to see how different data processing techniques affects accuracy. But those experiments do not belong in the final report, just include the final processing method and state if you did some experimentation to arrive at that method.

1. Train machine learning models to predict drug response on each cell line dataset (GDSC_Sanger, GDSC_Broad, CCLE) using 10-fold cross-validation. You can use any machine learning model available to you in R and tune its performance. You would train a model for each drug and each cell line dataset. For those best, performaing models, make a table or heatmap (with values printed on the cells) for the CV R-squared values. columns: CCLE, GDSC_Broad, GDSC_Sanger and rows: ALPELISIB, BUPARLISIB, CETUXIMAB, ERLOTINIB, RIBOCICLIB, TRAMETINIB. Do not change this row/column order.

2. By now you have a tuned model for each drug and cell line dataset combination. What happens when you use a model trained on one cell line dataset on on the other two? For example, for a given drug say Cetuximab, apply the CCLE model on GDSC_Sanger and GDSC_Broad datasets and record R-squared for those predictions. Then, you can take the GDSC_Broad model and apply it on other two and record R-squared values. If you keep doing this, for each drug you will have six values for R-squared. Output median and mean values as a table, and make boxplots for R-squared values of each drug as obtained above.

3. Take the best model from sub-task 1 and also the combine the models models in an ensemble learner ( this could be simply averaging their output) and predict drug response on the PDX data set. Plot the metrics (R-squared, Pearson correlation and Spearman’s Correlation) for observed vs predicted efficacy of each drug. Did the simple model perform better than ensemble model?

• HINT : How to decide the best model from step-1 for a given drug, this could be simply looking at CV R-squared values, or pick the model which has similar CV R-squared values to the out-of-sample R-squared values on the two other datasets. Anything reasonable is OK.

• HINT : For the ensemble model, you can simply average output or re-train an ensembl model with cross-validation. You don’t have to use all the models in the ensemble, you can choose to remove a model if you believe it performs poorly. Anything is OK but not every approach will be performant. Use your imagination.

• IMPORTANT NOTE : the PDX data set is for test ONLY and it shouldn’t be used for training the model or optimizing any training parameter including ensemble models.

4. Compare the best performing models for PDX prediction from TASK 1 and TASK2. Make informative figures and briefly discuss (2-3 sentences) if multi-omics or panel sequencing is better for modeling response judging from the PDX results.

TASK 3: Improving Drug Response Predictions by Transforming the Data in a Batch-Aware Manner.

The PDX and the cell line datasets are biologically different. On top of this, technical differences due to data acquisition and processing exists. These differences might drive the patterns we observe more than the biology that is related to the drug response, and the differences create so-called “batches” of data where data points cluster by technical effects rather than patters of interesting biology. Ideally, we can try to remove these differences and train models on datasets where these differences are removed. In many cases, we may not be able to completely remove these effects or removing these effects could also remove interesting biological patterns. However, there is also potential to produce better models that could have higher performance on PDX models.

Keeping these in mind, the sub-tasks are as follows:

1. Apply batch-removal methods on multi-omics datasets (TASK 1) including PDX. You can choose to work with a subset or only one of the cell line datasets but PDX dataset must be in there prior to batch effect removal. After removal of the effect, train drug response prediction onlny on the cell line dataset(s) and predict on PDX. Report if R-squared improved over TASK1 or not. You can try multiple batch removal methods. If that’s the case, report the PDX prediction performance for each of the methods. Make easy to undetstand & presentable plots.

2. Apply batch-removal methods to panel sequnecing datasets (TASK 2) including PDX. You can choose to work with a subset or only one of the cell line datasets but PDX dataset must be in there prior to batch effect removal. After removal of the effect, train drug response prediction only on the cell line dataset(s) and predict on PDX. Report if R-squared improved over TASK 2 or not. You can try multiple batch removal methods. If that’s the case, report the PDX prediction performance for each of the methods. Make easy to undetstand & presentable plots.

• HINT : Here are a few things to try for batch effect removal. You don’t have to try these things, there are other methods but if you need a place to start these are possibly viable options.
  • use SVD, find eigenvectors associated with batch class, remove them and reconstruct data.
  • use linear regression to remove effects example: https://rdrr.io/github/LTLA/batchelor/man/rescaleBatches.html or https://rdrr.io/github/LTLA/batchelor/man/regressBatches.html
  • use single-cell RNA-seq methods that uses mutual-nearest neighbours to correct for batch effects: https://rdrr.io/github/LTLA/batchelor/man/fastMNN.html

You might have problems applying these methods to binarized datasets such as mutation matrices. But you can try data transformations to make them work. One easy trick to try on binary matrices is to set 0s to -1. Another one could be to try logistic regression residuals instead of linear regression residuals ( regressBatches function uses linear regression residuals). You should probably try to use these methods in each omic data type seperately and construct corrected matrices for each of them. There are methods that could work on multi-omics datasets to do matrix factorization, you can try them as well and use them in an SVD-like manner as suggested above, ex: find latent factors associated with batch and remove them.

Summary of Data to be Used in TASK 1

Below some of the necessary summary statistics about the data you will use to complete tasks. A summary of drug response data is given at the end of the document.

#read the data
CCLE_gex <- read.delim("data/CCLE/all_intersect_genes/gex.tsv") %>% 
  rename(cell_lines = X) #Change the name of cell lines column to cell_lines

#Check for the normality
 hist(unlist(CCLE_gex[,2:50]),xlab="Gene Expression values",
     main="The CCLE Gene Expression ",border="white",
     col="cornflowerblue",bins = 50)

#Plot the boxplots for the Gene Expression values of the first 50 columns
 boxplot(CCLE_gex[,2:50],outline=FALSE, col="cornflowerblue",
        xlab = "Gene Expression values",xaxt="n",
        ylab = "Frequency",main = "The CCLE Gene Expression")

#check the missing values of the data set
CCLE_gex[CCLE_gex == "NA"] <- NA
table(is.na(CCLE_gex)) #FALSE:23954580  
#There are no missing values

#read the data
CCLE_cnv <- read.delim("data/CCLE/all_intersect_genes/cnv.tsv") %>%
  rename(cell_lines = X) #Change the name of cell lines column to cell_lines

cols <- colorRampPalette(brewer.pal(6,name="Blues"))(12)

pheatmap::pheatmap(CCLE_cnv[,2:50],cluster_rows=FALSE,show_rownames = FALSE,main="Copy Number Data from the CCLE",color = cols)

#check the missing values of the data set
CCLE_cnv[CCLE_cnv == "NA"] <- NA
table(is.na(CCLE_cnv)) #FALSE:32031220   
#There are no missing values

CCLE_mutdelim <- read.delim(file = "data/CCLE/all_intersect_genes/mut.tsv", 
                               sep = "\t", header = TRUE, fill = TRUE,
                            comment.char = "#") #read the mutation data 

CCLE_mutdelim$Tumor_Seq_Allele2 <- CCLE_mutdelim$Tumor_Seq_Allele1
CCLE_mut <- read.maf(maf= CCLE_mutdelim) #to convert and read it as maf file

## -Validating
## --Removed 6264 duplicated variants
## --Non MAF specific values in Variant_Classification column:
##   5UTR
##   5Flank
##   Start_Codon_Del
##   Stop_Codon_Ins
##   3UTR
## -Silent variants: 355334 
## -Summarizing
## --Possible FLAGS among top ten genes:
##   TTN
##   MUC16
##   OBSCN
##   SYNE1
##   USH2A
## -Processing clinical data
## --Missing clinical data
## -Finished in 40.0s elapsed (39.5s cpu)

#Plot the summary of mutation data 
plotmafSummary(maf = CCLE_mut, rmOutlier = TRUE, addStat = 'median',
               dashboard = TRUE,titvRaw = FALSE)

#check the missing values of the data set
CCLE_mutdelim[CCLE_mutdelim == "NA"] <- NA
table(is.na(CCLE_mutdelim)) #FALSE: 43439347 , TRUE: 886853 
#There are missing values

#read the data
BROAD_gex <- read.delim("data/GDSC_broad/all_intersect_genes/gex.tsv") %>% 
  rename(cell_lines = X) #Change the name of cell lines column to cell_lines

#Check for the normality
 hist(unlist(BROAD_gex[,2:50]),xlab="Gene Expression values",
     main="The GDSC BROAD Histogram ",border="white",
     col="Green",bins = 50)

#Plot the boxplots for the Gene Expression values of the first 50 columns
 boxplot(BROAD_gex[,2:50],outline=FALSE, col="Green",
        xlab = "Gene Expression values",xaxt="n",
        ylab = "Frequency",main = "The GDSC BROAD BoxPlot")

#check the missing values of the data set
BROAD_gex[BROAD_gex == "NA"] <- NA
table(is.na(BROAD_gex)) #FALSE:11252228   
#There are no missing values

#read the data
BROAD_cnv <- read.delim("data/GDSC_broad/all_intersect_genes/cnv.tsv") %>%
  rename(cell_lines = X) #Change the name of cell lines column to cell_lines

cols <- colorRampPalette(brewer.pal(6,name="Greens"))(12)

pheatmap::pheatmap(BROAD_cnv[,2:50],cluster_rows=FALSE,show_rownames = FALSE,
                   main="Copy Number Data from the GDSC BROAD",col = cols)

#check the missing values of the data set
BROAD_cnv[BROAD_cnv == "NA"] <- NA
table(is.na(BROAD_cnv)) #FALSE:11252228   
#There are no missing values

BROAD_mutdelim <- data.table::fread(
  file = "data/GDSC_broad/all_intersect_genes/mut.tsv",sep = "\t", 
  header = TRUE, fill = TRUE) #read the mutation data

genesum <- maftools::getGeneSummary(MAF(BROAD_mutdelim)) %>% group_by(Hugo_Symbol) %>% select(Hugo_Symbol,AlteredSamples) %>%
  arrange() %>% head(10)

## -Processing clinical data
## --Missing clinical data

Variant_Type <- BROAD_mutdelim %>% group_by(Variant_Type) %>%
count() %>% rename(frequency = n)

var_plot <- ggplot(data=Variant_Type, aes(x=Variant_Type, y = frequency, fill=Variant_Type))+
  geom_bar(stat="identity") +
  coord_flip() +
  ggtitle("Variant Type") +
  ylab("Frequency") +
  xlab("Variant Type")

ggplot(data=genesum, aes(x=reorder(Hugo_Symbol, AlteredSamples), y = AlteredSamples, fill = AlteredSamples))+
  geom_bar(stat="identity") +
  coord_flip() +
  ggtitle("Top 10 Mutated Genes Based on Altered Samples") +
  ylab("Number of Altered Samples") +
  xlab("Gene Names")+
  theme(legend.position = "right")

ggplotly(var_plot)

#check the missing values of the data set
BROAD_mutdelim[BROAD_mutdelim == "NA"] <- NA
table(is.na(BROAD_mutdelim)) #FALSE: 22905359  , TRUE: 24115732  
#There are missing values

#read the data
SANGER_gex <- read.delim("data/GDSC_sanger/all_intersect_genes/gex.tsv") %>% 
  rename(cell_lines = X) #Change the name of cell lines column to cell_lines

#Check for the normality
 hist(unlist(SANGER_gex[,2:50]),xlab="Gene Expression values",
     main="The GDSC SANGER Histogram ",border="white",
     col="Orange",bins = 50)

#Plot the boxplots for the Gene Expression values of the first 50 columns
 boxplot(SANGER_gex[,2:50],outline=FALSE, col="Orange",
        xlab = "Gene Expression values",xaxt="n",
        ylab = "Frequency",main = "The GDSC SANGER BoxPlot")

#check the missing values of the data set
SANGER_gex[SANGER_gex == "NA"] <- NA
table(is.na(SANGER_gex)) #FALSE:7048704   
#There are no missing values

#read the data
SANGER_cnv <- read.delim("data/GDSC_sanger/all_intersect_genes/cnv.tsv") %>%
  rename(cell_lines = X) #Change the name of cell lines column to cell_lines
cols <- colorRampPalette(brewer.pal(6,name="Oranges"))(12)

pheatmap::pheatmap(SANGER_cnv[,2:50],cluster_rows=FALSE,show_rownames = FALSE,
                     main="Copy Number Data from the GDSC SANGER", col=cols)

#check the missing values of the data set
SANGER_cnv[SANGER_cnv == "NA"] <- NA
table(is.na(SANGER_cnv)) #FALSE:16355196    
#There are no missing values

SANGER_mutdelim <-data.table::fread(
  file = "data/GDSC_sanger/all_intersect_genes/mut.tsv",sep = "\t", 
  header = TRUE, fill = TRUE)#read the mutation data 

genesum <- maftools::getGeneSummary(MAF(SANGER_mutdelim)) %>%
  group_by(Hugo_Symbol) %>% select(Hugo_Symbol,AlteredSamples) %>% arrange() %>%
  head(10)

## -Processing clinical data
## --Missing clinical data

Variant_Type <- SANGER_mutdelim %>% group_by(Variant_Type) %>% count() %>%
  rename(frequency = n)

var_plot <- ggplot(data=Variant_Type, aes(x=Variant_Type, 
                                          y = frequency, fill=Variant_Type))+
  geom_bar(stat="identity") +
  coord_flip() +
  ggtitle("Variant Type") +
  ylab("Frequency") +
  xlab("Variant Type")

ggplot(data=genesum, aes(x=reorder(Hugo_Symbol, AlteredSamples), y = AlteredSamples, fill = AlteredSamples))+
  geom_bar(stat="identity") +
  coord_flip() +
  ggtitle("Top 10 Mutated Genes Based on Altered Samples") +
  ylab("Number of Altered Samples") +
  xlab("Gene Names")+
  theme(legend.position = "right")

ggplotly(var_plot)

#check the missing values of the data set
SANGER_mutdelim[SANGER_mutdelim == "NA"] <- NA
table(is.na(SANGER_mutdelim)) #FALSE: 24111028   , TRUE: 24115732   
#There are missing values

#read the data
PDX_gex <- read.delim("data/PDX_all/all_intersect_genes/gex.tsv") %>% 
  rename(cell_lines = X) #Change the name of cell lines column to cell_lines

#Check for the normality
 hist(unlist(PDX_gex[,2:50]),xlab="The PDX Dataset",
     main="The PDX Gene Expression ",border="white",
     col="Purple")

#Plot the boxplots for the Gene Expression values of the first 50 column
PDX_gex[,2:50] %>% boxplot(names = FALSE,outline=FALSE, col="Purple",
                           xlab = "Gene Expression values",ylim = c(0,400),
                           ylab = "Frequency",main = "The PDX Gene Expression")

#check the missing values of the data set
PDX_gex[PDX_gex == "NA"] <- NA
table(is.na(PDX_gex)) #FALSE:23954580   
#There are no missing values

#read the data
PDX_cnv <- read.delim("data/PDX_all/all_intersect_genes/cnv.tsv") %>%
  rename(cell_lines = X) #Change the name of cell lines column to cell_lines

cols <- colorRampPalette(brewer.pal(6,name="Purples"))(12)
pheatmap::pheatmap(PDX_cnv[,2:50],cluster_rows=FALSE,show_rownames = FALSE,
                     main="Copy Number Data from the PDX",color = cols)

#check the missing values of the data set
PDX_cnv[PDX_cnv == "NA"] <- NA
table(is.na(PDX_cnv)) #FALSE:16355196     
#There are no missing values

PDX_mutdelim <- read.delim(file = "data/PDX_all/all_intersect_genes/mut.tsv", 
                               sep = "\t", header = TRUE, fill = TRUE,
                            comment.char = "#") #read the mutation data 

hist(unlist(PDX_mutdelim[,84]),xlab="Mutation Values",
     main="The PDX Mutation Frequency",border="white",
     col="Purple",bins = 10)

#check the missing values of the data set
PDX_mutdelim[PDX_mutdelim == "NA"] <- NA
table(is.na(PDX_mutdelim)) #FALSE: 806164 , TRUE: 8432975    
#There are missing values

Summary of Data to be Used in TASK 2

Below are all the necessary summary statistics about the data you will use to complete tasks. A summary of drug response data is given at the end of the document.

#read the data
CCLE_cnv_pan <- read.delim("data/CCLE/oncoKB_intersect_genes/cnv.tsv") %>%
  rename(cell_lines = X) #Change the name of cell lines column to cell_lines

cols <- colorRampPalette(brewer.pal(6,name="Blues"))(12)

pheatmap::pheatmap(CCLE_cnv_pan[,2:50],cluster_rows=FALSE,show_rownames = FALSE,main="Copy Number Data from the CCLE",color = cols)

#check the missing values of the data set
CCLE_cnv_pan[CCLE_cnv_pan == "NA"] <- NA
table(is.na(CCLE_cnv_pan)) #FALSE:1540835    
#There are no missing values

CCLE_mutdelim_pan <- read.delim(file = "data/CCLE/oncoKB_intersect_genes/mut.tsv",
                               sep = "\t", header = TRUE, fill = TRUE,
                            comment.char = "#") #read the mutation data 

CCLE_mutdelim_pan$Tumor_Seq_Allele2 <- CCLE_mutdelim_pan$Tumor_Seq_Allele1
CCLE_mut_pan <- read.maf(maf= CCLE_mutdelim_pan) #to convert and read it as maf file

## -Validating
## --Removed 522 duplicated variants
## --Non MAF specific values in Variant_Classification column:
##   Stop_Codon_Ins
##   3UTR
##   Start_Codon_Del
##   5Flank
## -Silent variants: 31994 
## -Summarizing
## -Processing clinical data
## --Missing clinical data
## -Finished in 2.250s elapsed (2.000s cpu)

#Plot the summary of mutation data 
plotmafSummary(maf = CCLE_mut_pan, rmOutlier = TRUE, addStat = 'median',
               dashboard = TRUE,titvRaw = FALSE)

#check the missing values of the data set
CCLE_mutdelim_pan[CCLE_mutdelim_pan == "NA"] <- NA
table(is.na(CCLE_mutdelim_pan)) #FALSE: 3719170 , TRUE: 77150  

## 
##   FALSE    TRUE 
## 3719170   77150

#There are missing values

#read the data
BROAD_cnv_pan <- read.delim("data/GDSC_broad/oncoKB_intersect_genes/cnv.tsv") %>%
  rename(cell_lines = X) #Change the name of cell lines column to cell_lines

cols <- colorRampPalette(brewer.pal(6,name="Greens"))(12)

pheatmap::pheatmap(BROAD_cnv_pan[,2:50],cluster_rows=FALSE,
                   show_rownames = FALSE,
                   main="Copy Number Data from the GDSC BROAD",col = cols)

#check the missing values of the data set
BROAD_cnv_pan[BROAD_cnv_pan == "NA"] <- NA
table(is.na(BROAD_cnv_pan)) #FALSE:786753    
#There are no missing values

BROAD_mutdelim_pan <- data.table::fread(
  file = "data/GDSC_broad/oncoKB_intersect_genes/mut.tsv",sep = "\t", 
  header = TRUE, fill = TRUE) #read the mutation data

genesum <- maftools::getGeneSummary(MAF(BROAD_mutdelim_pan)) %>% group_by(Hugo_Symbol) %>% select(Hugo_Symbol,AlteredSamples) %>%
  arrange() %>% head(10)

## -Processing clinical data
## --Missing clinical data

Variant_Type <- BROAD_mutdelim_pan %>% group_by(Variant_Type) %>%
count() %>% rename(frequency = n)

var_plot <- ggplot(data=Variant_Type, aes(x=Variant_Type, y = frequency, fill=Variant_Type))+
  geom_bar(stat="identity") +
  coord_flip() +
  ggtitle("Variant Type") +
  ylab("Frequency") +
  xlab("Variant Type")

ggplot(data=genesum, aes(x=reorder(Hugo_Symbol, AlteredSamples), y = AlteredSamples, fill = AlteredSamples))+
  geom_bar(stat="identity") +
  coord_flip() +
  ggtitle("Top 10 Mutated Genes Based on Altered Samples") +
  ylab("Number of Altered Samples") +
  xlab("Gene Names")+
  theme(legend.position = "right")

 ggplotly(var_plot)

#check the missing values of the data set
BROAD_mutdelim_pan[BROAD_mutdelim_pan == "NA"] <- NA
table(is.na(BROAD_mutdelim_pan)) #FALSE: 2143308   , TRUE: 2256126  
#There are missing values

#read the data
SANGER_cnv_pan <- read.delim("data/GDSC_sanger/oncoKB_intersect_genes/cnv.tsv") %>%
  rename(cell_lines = X) #Change the name of cell lines column to cell_lines

cols <- colorRampPalette(brewer.pal(6,name="Oranges"))(12)

pheatmap::pheatmap(SANGER_cnv_pan[,2:50],cluster_rows=FALSE,show_rownames = FALSE,
                     main="Copy Number Data from the GDSC SANGER", col=cols)

#check the missing values of the data set
SANGER_cnv_pan[SANGER_cnv_pan == "NA"] <- NA
table(is.na(SANGER_cnv_pan)) #FALSE:786753     
#There are no missing values

SANGER_mutdelim_pan <-data.table::fread(
  file = "data/GDSC_sanger/oncoKB_intersect_genes/mut.tsv",sep = "\t", 
  header = TRUE, fill = TRUE)#read the mutation data 

genesum <- maftools::getGeneSummary(MAF(SANGER_mutdelim_pan)) %>%
  group_by(Hugo_Symbol) %>% select(Hugo_Symbol,AlteredSamples) %>% arrange() %>%
  head(10)

## -Processing clinical data
## --Missing clinical data

Variant_Type <- SANGER_mutdelim_pan %>% group_by(Variant_Type) %>% count() %>%
  rename(frequency = n)

var_plot <- ggplot(data=Variant_Type, aes(x=Variant_Type, 
                                          y = frequency, fill=Variant_Type))+
  geom_bar(stat="identity") +
  coord_flip() +
  ggtitle("Variant Type") +
  ylab("Frequency") +
  xlab("Variant Type")

ggplot(data=genesum, aes(x=reorder(Hugo_Symbol, AlteredSamples), y = AlteredSamples, fill = AlteredSamples))+
  geom_bar(stat="identity") +
  coord_flip() +
  ggtitle("Top 10 Mutated Genes Based on Altered Samples") +
  ylab("Number of Altered Samples") +
  xlab("Gene Names")+
  theme(legend.position = "right")

 ggplotly(var_plot)

#check the missing values of the data set
SANGER_mutdelim_pan[SANGER_mutdelim_pan == "NA"] <- NA
table(is.na(SANGER_mutdelim_pan)) #FALSE: 2143308 , TRUE: 2256126    
#There are missing values

#read the data
PDX_cnv_pan <- read.delim("data/PDX_all/oncoKB_intersect_genes/cnv.tsv") %>%
  rename(cell_lines = X) #Change the name of cell lines column to cell_lines

cols <- colorRampPalette(brewer.pal(6,name="Purples"))(12)
pheatmap::pheatmap(PDX_cnv_pan[,2:50],cluster_rows=FALSE,show_rownames = FALSE,
                     main="Copy Number Data from the PDX",color = cols)

#check the missing values of the data set
PDX_cnv_pan[PDX_cnv_pan == "NA"] <- NA
table(is.na(PDX_cnv_pan)) #FALSE:352317      
#There are no missing values

PDX_mutdelim_pan <- read.delim(file = "data/PDX_all/oncoKB_intersect_genes/mut.tsv", 
                               sep = "\t", header = TRUE, fill = TRUE,
                            comment.char = "#") #read the mutation data 

hist(unlist(PDX_mutdelim_pan[,84]),xlab="Mutation Values",
     main="The PDX Mutation Frequency",border="white",
     col="Purple",bins = 10)

#check the missing values of the data set
PDX_mutdelim_pan[PDX_mutdelim_pan == "NA"] <- NA
table(is.na(PDX_mutdelim_pan)) #FALSE: 75781 , TRUE: 787536     
#There are missing values